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Abstract 

The thermodynamical stability of a set of circular double helical molecules is analyzed by path 
integral techniques. The minicircles differ only m. i) the radius and ii) the number of base pairs (iV) 
arranged along the molecule axis. Instead, the rise distance is kept constant. For any molecule size, 
the computational method simulates a broad ensemble of possible helicoidal configurations while 
the partition function is a sum over the path trajectories describing the base pair fluctuational 
states. The stablest helical repeat of every minicircle is determined by free energy minimization. 
We hnd that, for molecules with N larger than 100, the helical repeat grows linearly with the 
size and the twist number is constant. On the other hand, by reducing the size below 100 base 
pairs, the double helices sharply unwind and the twist number drops to one for N = 20. This is 
predicted as the minimum size for the existence of helicoidal molecules in the closed form. The 
helix unwinding appears as a strategy to release the bending stress associated to the circularization 
of the molecules. 

PACS numbers: 87.14.gk, 87.15.A-, 87.15.Zg, 05.10.-a 
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Introduction 


A large amount of physical and biochemical research on circular DNA has been carried 
out over the last hfty years following the hrst observations that the single stranded DNA 


of coliphage 0X174 was a ring 
closed circular duplex molecules 


and the tumor virus DNA was in the form of covalently 


Being less degradable than linear DNA, circular DNA 
has recently revealed a potential for those biotechnological applications such as assembly of 
nanotubes which require both rigidity and flexibility of the building blocks QQ. 

Bending flexibility of short DNA fragments has been appreciated since it was pointed 
out jo] that linear molecules of order 100 base pairs (bps) can spontaneously close into 


ler cyclization probabilities than 
7H9j. The circularization of a DNA 


circles, even in the absence of proteins, with much hig 
those predicted by conventional worm like chain models 
fragment causes a bending stress in the loop that may be released by a local untwisting of 
the double helix with formation of fluctuational bubbles comprising a few transiently broken 
bps While bubbles are essential to key biological functions such as replication and 

transcription of the genetic code 


llMlSl]. their length is related to the degree of supercoiling 
as shown by atomic force microscopy imaging of circular plasmids Q. 

In general, bending and twisting degrees of freedom are intertwined in DNA and their 
interplay crucially depends on the molecule size. 

Recently, short sequences of extra-chromosomal circular DNA have been found, both in 
mouse tissues and human cells, in a range of 80 — 2000 bps with length distributions peaked 
between 200 and 400 bps 17|]. These hndings lend support to the prediction that there 
should be a minimum size for the existence of double stranded loops whereas, below that 
size, circular molecules can be only single stranded 18| . However, the determination of such 


threshold is still an open question also in view of the fact that duplex DNA minicircles have 


been obtained, by a modihed ligase-assisted protocol 


19| . even in the 60 bps range. 


Specihcally in the latter, disruptions of the regular helical structure have been detected by 
single-strand-specihc endonuclease experiments carried out on DNA minicircles with various 


radii 


lOj. Such local disruptions may appear either in the form of kinks which move one 


base pair out of the stack [20|, l2]| or as a breaking of a few base pair hydrogen bonds j22| . 
the latter event being energetically more costly. 

These issues have been attacked by a theoretical model based on the path integral for- 
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malism 23|] and applied to three minicircles having comparable relative contents of AT- and 


GC-bps, albeit with different radii 


24j |: analysis of the bnbble prohles has shown that base 


pair disruptions and helix unwinding are in fact more pronounced in the smallest minicircle 
with 66 bps (corresponding to a bending angle of ~ 6°) while larger loops are relatively 
stabler. Whereas these hndings conhrm the soundness of the path integral method, it is 
still unsettled whether a critical radius of curvature for the appearance of bps disruptions in 
small circles indeed exists and to which extent such critical value may depend on the DNA 
sequence specihcities. 

To shed light on these issues we investigate in this work the thermodynamical stability 
of a broader set of molecules spanning the range between 20 and 140 bps. In view of the 
shortness of these systems we assume that, for any given length, the cyclization occurs 


with hnite probability and the molecules are thus in the circular form 


25| . Furthermore, 


we admit that the minicircles could be found either in the single or in the double stranded 
helicoidal conhguration, keeping the model general enough to deal with both cases. For each 
minicircle size, we simulate a large ensemble of possible twist conhgurations and compute 
the free energy prohles as a function of the average helical repeat, i.e. the number of bps per 
helix turn. This methods permits to determine the energetically most convenient helicoidal 
conformation for a molecule that, by virtue of the circularization process, has incorporated 
a certain amount of bending stress. 


Model 


In general, the quantum statistical partition function of a system is determined by per¬ 
forming an analytic continuation of the quantum mechanical partition function to the imag¬ 


inary time r ^ U, t be.ng the tea. t.me M- r .s defined with.n a ^nge whose upper bound 

is given by the inverse temperature as in the Matsubara method [27|. Accordingly, the Eu¬ 
clidean action takes the place of the mechanical canonical action and, in the path integral 
formalism 28|, the statistical partition function is an integral in the conhguration space 


over (closed) paths running along an imaginary time axis. The largest contribution to the 
statistical partition function comes from those trajectories for which the sum of the kinetic 
and potential energy is small. The classical partition function is obtained by replacing the 
quantum thermal wavelength with the classical one as described in detail in 29 |. 
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The focus of our investigation is a system made of a homogeneous set of N free purine- 
pyrimidine bps with reduced mass /i. Let’s dehne r* as the inter-strand fluctuational vector 
describing the dynamics of the stretching mode for the i — th base pair. The application of 
the path integral method to our classical system relies on the idea that, at hnite temperature, 
the base pair displacements can be thought of as trajectories, Xi^r): 


± |ri| Xi{T); re [0 ,/?] , 


( 1 ) 


where (3 = (ksT) ks is the Boltzmann constant and T is the temperature. As the 
trajectories are closed, a;i(0) = Xi{(3), they can be expanded in Fourier series: 


Xi{T) = (xo)* ^ cos(a;mr) {b^)iSm{u^T) 

m=l 

2mn 


The set of coefficients {(xq)*, {am)i, represents a point in the path conhguration 

space corresponding to a ffuctuational state for the i — th base pair. Consistently one defines, 
over the space of the Fourier coefficients, an integration measure f Dxi which normalizes 
the partition function Zk for the freely fluctuating N particles: 


^ I 

Zk= Wji Dxi exp [-A, 

eJ “r 

i=l >^0 




'a:l = 1 


Ak\x] = 


k' •2r \ 


(3) 


Ac is the classical thermal wavelength and Ak[x\ is the kinetic action which, as a conse¬ 
quence of the real space mapping, is a dimensionless quantity. Hence, Eq. (|3]) sets the zero 
for the free energy, InZ^, of the system. 

Next, as shown in Fig. [H we uniformly arrange the N objects (represented by the blue 
points) on a circle, with radius i?, which represents the bent molecule axis. The i — th object 
can fluctuate with respect to R and its vector (r,) describes a (red shaded) orbit lying on 
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a plane which is bent by with respect to the (x, y) plane. In the local reference system, 
centered in Oj, the i — th object makes a positive twist angle (Ot) with respect to the x' 
axis and adjacent objects along the stack are twisted by a constant angle. The rise distance 
between neighbor objects along the circnmference is: d = 2tiR/N . Thns we have bnilt a 
circnlar model for a set of particles arranged in a helicoidal conformation. The model can 
be applied to both a single stranded and a donble stranded molecnle. In the former case, r* 
represents the i — th base along the strand. In the latter, r* is the inter-strand flnctnational 
vector, for the i — th base pair, with respect to the eqnilibrinm position which corresponds 
to the minimnm for the hydrogen bond potential (see below). In the following we will refer 
to a set of N bps with the caveat that very short circnlar molecnles may more likely be 
fonnd in the single strand conhgnration. The valnes d = 3.4A and /i = 300 amn pecnliar 
of DNA molecnles are hereafter taken. 

Accordingly, with respect to the O reference system (see Fig. [1]), the general vector tj for 
the i — th base pair is: 


[ti)^= |ri| cos^icos^i 
{ti)y = (i?+|rj| sin^j) cos0i 

{ti)^ = {R+ |ri|)sin0j. (4) 


The x-axis is normal to the sheet plane and R lies on the (y, z) plane. For small size 
molecnles, it can be reasonably assnmed that the relative bending between adjacent orbital 
planes along the stack is constant and depends only on the circle length. Then, the bending 
angle is, 0* = {i—l)27i/N+(j)s, while the twisting is measnred by, 9i = {i—l)27i/h+9s, where 
h is the nnmber of bps per helix tnrn [^, M]- The compntation snms over a distribntion of 
(f)s and 9s in order to avoid to pin the hrst base pair in the seqnence to a specihc angnlar 
position. Fnrthermore, as the molecnle axis lies on a plane, the Writhe (measnring the 


spatial coiling of the axis itself 


32[ | ) vanishes and the integer Lin 


closed-circnlar molecnle coincides with the Twist nnmber Tw 


c ing nnmber {Lk) for the 
34|- The latter is given 
hj Tw = N/h which, however, is generally not an integer. While in long molecnles one 
can take Lk as the closest integer to N/h, the extra twist reqnired to join the strand ends 
becomes large in short seqnences thus signihcantly reducing the cyclization probability j?]. 
Accordingly, as h is here a variable, we assume to deal with closed-circular molecules only if 
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FIG. 1: (Color online) Schematic of the base pairs (blue dots) stacked along the circular molecule 
backbone whose radius R lies on the {y, z) plane, is the inter-strand fluctuational vector, for the 
i — th base pair, spanning the red shaded orbit. For |rj| = 0, the orbit shrinks into the fluctuations 
free state. Adjacent orbits along the molecule axis are bent by an angle (j)i. 9i is the twist angle 
for the i — th base pair as shown in the reference system Oj. 

N/h is an integer and, in these cases, also closure conditions on the radial fluctuations are 
implemented in the computation, |rAr+i| = |ri|. Whenever such conditions are not fulhlled, 
our circular molecules are open. 

As the efficacy of the method crucially depends on the partition of the Tw range, we take 
a small incremental step, ATw = 0.0125, so that any molecule may assume in principle a 
large ensemble of possible twisted conformations, each with a specihc h. The computational 
task consists in evaluating the free energy prohles within a broad range of h values. It follows 
that the CPU time grows fast with the size of the molecule. 

In analogy with Eq. ([T]), the path integral method is applied to the bent conhguration. 
Accordingly, tj maps onto a time dependent path fluctuation: 


tj| ^ r]i{T) , 


(5) 
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which, in the circular geometry, reads: 


T]i{T) = [R^ + XiiT)^^ + 2R\Xi{T)\f{ei,(pi)Y^‘^ 

f{6i, (pi) = sin Oi cos^ (pi + sin^ (pi. (6) 

Then, the partition function for the ensemble of N bps oscillating around the fluctuations 
free circle is given by: 


N 

^0 = n f Dxiexp[-Ao[r]]] 

*=1 9 s As 


N 


Mv] 




i=l 


f/3 


^okl = 


dr 


A'' -2/ \ 
2 "■ 


(7) 


The action Ao[t]] accounts for the bending energy cost due to the deformation which 
closes the array of N free particles into a loop. This is proved by the following analytical 
argument. 

Take for simplicity only one Fourier component in Eq. (|2]) (and set {bi)i = 0). Then the 
path amplitude can be approximated by: Xi^r) ~ (a;o)i + (ai)j cos(a;ir) and Ao[rii\ in Eq. ([7]) 
reads: 


dr 


d' •2f \ 

2 (^) 


/S 




dr sin 




9i[r) = 


[|a;i(r)| + Rf{9i,(pi)]' 


[R^ + a:i(r)2 + 2i?|xi(r)|/(6'i, (pi)] 


isi)t ^ 


TT^{ai)l 


>^1 


( 8 ) 


As the action is not quadratic in the Fourier coefficients, Zq does not decouple into 
a product of Gaussian integrals. This marks the difference with respect to Zk in Eq. ([3|). 
Setting R = 0, the Gaussian approximation is recovered for any Fourier component, Ao[77i] ~ 
(si)^ and Zq —)■ Z^. 

In general, for a circular molecule with finite R, the normalization condition in Eq. ([3]) 
is not fulfilled. Hence, the free energy level for the loop, Fq = /9“^lnZo, is not zero. Fq is 
computed via Eq. ([7]) with the integration measure given in the first of Eqs. fl3l) . 
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A mesoscopic model for double stranded DNA molecules should incorporate in the action 
also the terms describing the effective hydrogen bond base pair interactions. This is done 
by including: a) the Morse potential, VM[^i('r)] = Di[exp{—bi{rii{r) — R)) — which 
also accounts for the repulsive electrostatic interaction between phosphate groups. Di is 
the pair dissociation energy while the inverse length bi sets the potential range. h) A 
solvent potential which enhances the barrier for base pair breaking and stabilizes the strands: 
Vsoi[Tlii.T)] = -A/s(tanh((? 7 i(r)-i?)//J-l) [3^. The factoris related to the counterions 
concentration in the solvent and Is tunes the width of the solvent barrier. Both terms and 
parameters choice, A = 45meld, b^ = 2 A“^, fg = 0.1, Ig = 0.5 A, are described in detail 
in 3^. Hence the partition function for our homogeneous ds-circular molecules reads: 


N 

Zi = JJ ^ (b Dxiexp[-Ai[r]]] 

i=l 6s,(l>s 


= 1 ds,(t>s 

N ./3 


= E / 

i.l ■'0 


dr 


^vK'r) + Vi[r]iiT)] 


= VM[Viir)] + VsoihiT)] ■ 


(9) 


and the free energy is computed via A = /d ^In A [36 1. 


Results and Discussion 


Given a molecule with N bps, for any twist conhguration, the programme sums over 
the sets of Fourier coefficients corresponding to ~ 10^ ffuctuational states for each base 
pair. This suffices to get numerical convergence in the partition function. The computation 
includes also large (of the order of a few Angstroms) path amplitudes Xi{T) which measure 
the distance between the pair mates on complementary strands. However, such amplitudes 
should not exceed the radius R in order to preserve the overall circular shape of the molecule. 
Technically, too large ffuctuational amplitudes can be discarded as they would yield a high 
kinetic action (see Eq. ([H])) hence, a small contribution to the partition function. 

The free energy prohles are shown in Figs. [2] for three circular molecules. A and A are 
compared in each panel. The calculations are carried out at room temperature being aware 
that thermal effects would loosen the base pair bonds, favor the molecules unwinding, ulti¬ 


mately leading to denaturation in the high temperature regime 


37| |. The input parameters 
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FIG. 2: (Color online) Room temperature free energies (in meV) for three helicoidal circular 
molecules as function of the helical repeat. Fq and Fi are computed by Eqs. d?]) and Eqs. ([9]), 
respectively. 
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of the model, i.e. rise distance, reduced mass and hydrogen bond potential, are common to 
all circles. 

While oscillating patterns are found in all the free energy versus h plots, the location 
of the free energy minimum strongly varies with the molecule size. This indicates that, for 
each minicircle, the thermodynamically stablest conhguration is associated to a specihc twist 
conformation. In particular, the smallest circle (among the three ones displayed) takes the 
equilibrium state with the largest helical repeat, h ~ 16, suggesting that the helix unwinding 
is more pronounced in small size molecules where the bending of the backbone is strong. 
Importantly, in each panel, the minima for Fi coincide with those for Fq signaling that the 
molecule unwinding is driven by the bending rather than by the hydrogen bonds, at least 
at room temperature. Then, the unwinding patterns of single and double stranded circular 
molecules are similar. On the other hand, the potential in Eqs. (jO]) contributes to order the 
molecular structure and, accordingly, the Fi values are shifted upwards on the energy scale 
with respect to Fq. 

Fig. |3] summarizes the results of the free energy computation for a set of circles with 
size N G [20, 140]. The helical repeat values which minimize the free energies are reported. 
Interestingly the curve is not monotonous and the most twisted conformation {h ~ 5) is 
found for N = 100. For 100 < N < 140, the h equilibrium values increase smoothly and 
linearly with N so that the twist number remains essentially constant, Tw = 20. This 
is however a property which should be investigated also for larger circles than those here 
examined. 

Instead, for N < 100, h sharply grows. This is physically consistent with the fact 
that short molecules, once bent, have small angles between adjacent bps along the stack. 
As a consequence, short circular molecules tend to untwist in order to attain an energetic 
conhguration which may efficiently release the bending stress. Furthermore, N/h is generally 
not an integer indicating that, in the very small size range, circular molecules are unlikely 
to be found in the closed conformation. 

Only for the loop with N = 20, one single turn of the helix involves all 20 bps that is, 
Tw = 1 and the bending angle is 18°. It follows that, in our model, no closed helicoidal 
molecules can be found under such size: a prediction which would be interesting to check 
experimentally. 

The helical repeats selected in Fig. [3] represent the most probable twist conformations 
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FIG. 3: (Color online) Helical repeat values obtained by minimization of the free energy of mini¬ 
circles with different size. 


associated to the deepest minimum in the free energy prohles with multiple wells. It should 
be however remarked that thermal effects may transiently drive every molecule out of its 
deepest well bringing it to one of the energetically close states associated to other h values 
38|. 


In conclusion, the picture emerging from this path integral computation points to the 
existence of a critical size of order 100 bps below which the molecule helix progressively 
unwinds. The applied model is based on the actions in Eq. ([7]) and Eq. ([2]) accounting for the 
kinetic energy and the stabilizing hydrogen bond potential. A more structured description 
including also the intra-strand base pair stacking may certainly rehne the predictions for the 
most probable helical repeat values of specihc sequences but it is not expected to change 


the trend of the results for the helical unwinding versus N presented in this work 




Thus, applying a criterion of thermodynamical stability, we have found that the unwind¬ 
ing of small helicoidal circular molecules hosting a uniform distribution of nucleotides is 
driven by the size of the loop and it is not essentially affected by the sequence specihcities. 
The helix unwinds in order to release the stress associated to the bending of the molecule 
backbone. The equilibrium twist conformations have been determined by minimizing the 
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molecules free energy under the assumption that the bending angle between adjacent base 
pairs along the stack is constant. For circles larger than those here examined, such as¬ 
sumption may not hold and, accordingly, bending fluctuational effects 


40 


41| should be 


incorporated in the theoretical description. Finally, for short double stranded molecules, 
a direct computation of the cyclization probabilities in the framework of the path integral 
method may contribute to clarify whether a threshold size for closed circular conformations 
indeed exists. 
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